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Abstract 

In a diffusion-based molecular communication network, transmitters and receivers com¬ 
municate by using signalling molecules (or ligands) in a fluid medium. This paper assumes 
that the transmitter uses different chemical reactions to generate different emission patterns 
of signalling molecules to represent different transmission symbols, and the receiver consists 
of receptors. When the signalling molecules arrive at the receiver, they may react with the 
receptors to form ligand-receptor complexes. Our goal is to study the demodulation in this 
setup assuming that the transmitter and receiver are synchronised. We derive an optimal 
demodulator using the continuous history of the number of complexes at the receiver as the 
input to the demodulator. We do that by first deriving a communication model which in¬ 
cludes the chemical reactions in the transmitter, diffusion in the transmission medium and 
the ligand-receptor process in the receiver. This model, which takes the form of a continuous¬ 
time Markov process, captures the noise in the receiver signal due to the stochastic nature of 
chemical reactions and diffusion. We then adopt a maximum a posteriori framework and use 
Bayesian Altering to derive the optimal demodulator. We use numerical examples to illustrate 
the properties of this optimal demodulator. 

Keywords: Molecular communication networks; modulation; demodulation; maximum a posteri¬ 
ori; optimal detection; stochastic models; Bayesian filtering; molecular receivers. 
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1 Introduction 


Molecular communication is a promising approach to realise communications among nano-scale 
devices [DEIEIII]. There are many possible applications with these networks of nano-devices, for 
example, in-body sensor networks for health monitoring and therapy |S1 E]- This paper considers 
diffusion-based molecular communication networks. 

In a diffusion-based molecular communication network, transmitters and receivers communicate 
by using signalling molecules or ligands. The transmitter uses different time-varying functions 
of concentration of signalling molecules (or emission patterns) to represent different transmission 
symbols. The signalling molecules diffuse freely in the medium. When signalling molecules reach 
the receiver, they react with chemical species in the receiver to produce output molecules. The 
counts of output molecules over time is the receiver output signal which the receiver uses to decode 
the transmitted symbols. 

Two components in diffusion-based molecular communication system are modulation and de¬ 
modulation. A number of different modulation schemes have been considered in the literature. 
For example, imi consider Concentration Shift Keying (CSK) where different concentrations of 
signalling molecules are used by the transmitter to represent different transmission symbols. Other 
modulation techniques that have been proposed include Molecule Shift Keying (MSK) [H |9] , Pulse 
Position Modulation (PPM) [TU], Amplitude Shift Keying (ASK)[TT], Frequency Shift Keying (FSK) 
[T^ . and token communication [T3]. This paper assumes that the transmitter uses different chem¬ 
ical reactions to generate the emission patterns of different transmission symbols. The motivation 
to use this type of modulation mechanism is that chemical reactions are a natural way to produce 
signalling molecules, e.g. the papers ong study a number of molecular circuits (which are sets 
of chemical reactions) that can produce oscillating signals, and the paper [16] discusses a number 
of signalling mechanisms in living cells. 

We assume the receiver consists of receptors. When the signalling molecules (ligands) reach 
the receiver, they can react with the receptors to form ligand-receptor complexes (which are the 
output molecules in this paper). We consider the problem of using the continuous-time history of the 
number of complexes for demodulation assuming that the transmitter and receiver are synchronised. 
The ligand-receptor complex signal is a stochastic process with three sources of noise because the 
chemical reactions at the transmitter, the diffusion of signalling molecules and the ligand-receptor 
binding process are all stochastic. We derive a continuous-time Markov process (CTMP) which 
models the chemical reactions at the transmitter, the diffusion in the medium and the ligand- 
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receptor binding process. By using this model and the theory of Bayesian hltering, we derive the 
maximum a posteriori (MAP) demodulator using the continuous-time history of the number of 
complexes as the input. 

This paper makes two key contributions: (1) We propose to use a CTMP to model a molecular 
communication network with chemical reactions at the transmitter, a diffusive propagation prop¬ 
agation medium and receptors at the receiver. The CTMP captures all three sources of noise in 
the communication network. (2) We derive a closed-form expression for the MAP demodulation 
hlter using the proposed CTMP. The closed-form expression gives insight into the important ele¬ 
ments needed for optimal demodulation, these are the timings at which the receptor bindings occur, 
the number of unbound receptors and the mean concentration of signalling molecules around the 
receptors. 

The rest of the paper is organised as follows. Section [^discusses related work. Section [^presents 
the system assumptions, as well as a mathematical model from the transmitter to the ligand-receptor 
complex signal based on CTMP. We derive the MAP demodulator in Section and illustrate its 
numerical properties in Section Finally, Section concludes the paper. 

2 Related work 

There is a growing interest to understand molecular communication from the communication engi¬ 
neering point of view. For recent surveys of the held, see [D El El SI El . We divide the discussion 
under these headings: transmitters, receivers, models and others. 

Transmitters. A number of different types of transmission signals have been considered in the 
molecular communication literature. The papers HHIC] assume that the transmitter releases the 
signalling molecules in a burst which can be modelled as either an impulse or a pulse with a hnite 
duration. A recent work in [12] assumes that the transmitter releases the molecules according to a 
Poisson process. In this paper, we instead assume that the transmitter uses different sets of chemical 
reactions to generate different transmission symbols and we use CTMP to model these transmission 
symbols. Since a Poisson process can also be modelled by a CTMP, the transmission process in this 
paper is more general than that of [12]. Our CTMP model can also deal with an impulsive input 
by using an appropriate initial condition for the CTMP. The use of CTMP as an end-to-end model 
— which includes the transmitter, the medium and the receiver — does not appear to have been 
used before. 
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Receivers. Demodulation methods for diffusion-based molecular communication have been 
studied in [201 ED. Both papers also use the MAP framework with discrete-time samples of the 
number of output molecules as the input to the demodulator. Instead, in this paper, we consider 
demodulation using continuous-time history of the number of complexes. 

The demodulation from ligand-receptor signal has also been considered in [TS]. The key differ¬ 
ence is that [iHj uses a linear approximation of the ligand-receptor process while we use a non-linear 
reaction rate. 

The capacity of molecular communications based on ligand-receptor binding has been studied 
in [221 [25] assuming discrete samples of the number of complexes are available. A recent work 
considers the capacity of such systems in the continuous-time limit. Instead of focusing on the 
capacity, our work focuses on demodulation. 

Receiver design is an important topic in molecular communication and has been studied in many 
papers, some examples are [SSlEQlEniEIlEZ]. These papers either use one sample or a number of 
discrete samples on the count of a specihc molecule to compute the likelihood of observing a certain 
input symbols. This paper takes a different approach and uses continuous-time signals. 

Another approach of receiver design for molecular communication is to derive molecular circuits 
that can be used for decoding. An attempt is made in [T2] to design a molecular circuit that can 
decode frequency-modulated signals. However, the work does not take diffusion and reaction noise 
into consideration. A recent work in [28] analyses end-to-end molecular communication biological 
circuits from linear time-invariant system point of view. The work in [22] compares the information 
theoretic capacity of a number of different types of linear molecular circuits. This paper differs from 
the previous work in that it uses a non-linear ligand-receptor binding model. 

The noise property of ligand-receptor for molecular communication has been characterised in 
[50] . The case for non-linear ligand-receptor binding does not appear to have an analytical solution 
and EDI derives an approximate characterisation using a linear reaction rate assuming that the 
number of signalling molecules around the receptor is large. This paper uses a non-linear ligand- 
receptor binding model and no approximation is used in solving the hltering problem. 

Models. This paper uses the Reaction Diffusion Master Equation (RDME) [5T] framework to 
model the reactions and diffusion in the molecular communication networks. RDME assumes that 
time is continuous while the diffusion medium is discretised into voxels. This results in a CTMP 
with hnite number of (discrete) states. RDME has been used to model stochastic dynamics of cells 
in the biology literature [52]. An attraction of RDME is that it has the Markov property which 
means that one can leverage the rich theory behind Markov process. 
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The author of this paper has previously used an extension of the RDME model, called the 
RDME with exogenous input (RDMEX) model, to study molecular communication networks in 
[33l [3T[ l35l [36]. The RDMEX assumes that the times at which the transmitter emits signalling 
molecules are deterministic. This results in a stochastic process which is piecewise Markov or the 
Markov property only holds in between two consecutive emissions by the transmitter. In this paper, 
we assume the transmitter uses chemical reactions to generate the signalling molecules. Therefore, 
the emission timings are not deterministic but are governed by a stochastic process. 

In this paper, we assume that the propagation medium is discretised in the voxels. An alternative 
modelling paradigm that has been used in a number of molecular communication network papers 
[acHicn] is that the transmitter or receiver has a non-zero spatial dimension (commonly modelled 
by a sphere) while the propagation medium is assumed to be continuous. (Note that though |T9] 
does not explicitly state the dimension of the receiver, one can infer from the fact that the receiver 
must have a non-zero dimension because it has a non-zero probability of receiving the signalling 
molecules.) We believe the technique in this paper can be adapted to this alternative modelling 
paradigm and we do not expect this alternative modelling paradigm will change the results in this 


paper; we will explain this in Section 4.4.2 


There is a rich literature in the modelling of biological systems discussing the difference between: 
(1) The particle approach which has a continuous state space because the state of a particle is its 
position; and (2) The mesoscopic approach (the approach in this paper) which discretises the 
medium into discrete voxels and consider the number of molecules in the voxels as the state. The 
first approach is more accurate but the computation burden can be high 1371, while the second 
approach is accurate for appropriate discretisation [381 EH- There are also hybrid approaches too. 
An overview of various modelling and simulation approaches can be found in p7] . 

Others: The results of this paper may also be of interest to biologists who are interested to 
understand how living cells can distinguish between different concentration levels [39l HQ]. The 
result of this paper can be viewed as a generalisation of [39] which studies how cells can distinguish 
between two constant levels of ligand concentration. 


3 End-to-end communication models 

This paper considers diffusion-based molecular communication with one transmitter and one receiver 
in a fluid medium. Figure [T] gives an overview of the setup considered in this paper. The transmitter 
uses different chemical reactions to generate the emission patterns of different transmission symbols. 
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The transmitter acts as the source and emitter of signalling molecules. The signalling molecules 
diffuses in the fluid medium. The front-end of the receiver consists of a ligand-receptor binding 
process and the back-end consists of the demodulator with the number of complexes as its input. 


In this section, we hrst describe the system assumptions in Section 3.1, We then present, in 


Section 3.2, an end-to-end model which includes the transmitter, the transmission medium and the 
ligand-receptor binding process in the receiver, see the dashed box in Figure [T} The end-to-end 
model is a CTMP which includes chemical reactions in the transmitter, diffusion in the medium 
and the ligand-receptor binding process in the receiver. 


3.1 Model assumptions 

We assume that the medium (or space) is discretised into voxels while time is continuous. This 
modelling framework results in a RDME [3IISIII12], which is a CTMP commonly used to model 
systems with both diffusion and reactions. In addition, we assume the communication uses only 
one type of signalling molecule (or ligand) denoted by S. We divide the description of our model 
into three parts: transmission medium, transmitter and receiver. We begin with the transmission 
medium. Table summaries the frequently used notation and chemical symbols. 

3.1.1 Transmission medinm 

We model the transmission medium as a three dimensional (3-D) space and partition the space into 
cubic voxels of volume W^. Figure shows an example of a medium which has a dimension of 4 
voxels along both the x and ^/-directions, and 1 voxel in the 2 ;-direction. (Note that Figure [^should 
be viewed as a projection onto the x — y plane.) In general, we assume the medium to have N^, Ny 
and Nz voxels in the x, y and z directions where iVa,, Ny and W are positive integers. In Figure 
Nx = Ny = A and Nz = 1. We also use W = NxNyNz to denote the total number of voxels. 

We refer to a voxel by a triple {x, y, z) where x, y and are integers or by a single index 
^ e [1, Ny]. Figure shows the triples for some of the voxels. The single index ^ is calculated from 
the triple {x, y, z) by using ^{x, y,z) = x + Nx{y — 1) + NxNy[z — 1). The single indices for voxels 
are shown in the top right-hand corner of the voxels in Figure 

Diffusion is modelled by molecules moving from one voxel to a neighbouring voxel. For examples, 
in Figure molecules can diffuse from Voxel 1 to Voxels 2 or 5, from Voxel 2 to Voxels 1, 3 and 
6, and so on. The diffusion of molecules between neighbouring voxels is indicated by the two-way 
arrows in Figure]^ 
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Symbol 

Meaning 

W 

Dimension of one side of a voxel 

D 

Diffusion constant 

d 

Diffusion rate between neighbouring voxels 

N, 

Total number of voxels 

M 

Total number of receptors 

A 

Reaction rate constant for the binding reaction 

A 

\ _ A 

h 

Reaction rate constant for the unbinding reaction 

s 

A transmission symbol 

b{t) 

Number of complexes at time t 

riiit) 

Number of signalling molecules in voxel i at time t 

Nit) 

Equation Q. A vector containing the number of signalling molecules in each voxel, the 
counts of intermediate chemical species in the transmitter and the cumulative count of 

the number of molecules that have left the system 

Csit) 

The mean number of signalling molecules in the receiver voxel at time t if the trans¬ 
mitter sends symbol s 

U{t) 

The cumulative number of times the receptors have switched from the unbound to 


bound state at time t 

E 

An unbound receptor 

S 

A signalling molecule 

C 

A complex 


Table 1: Notation and chemical symbols 
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We assume that the signalling molecule S is the only diffusible chemical species in our model and 
the diffusion coefficient for S' is D. This means the signalling molecules diffuse from one voxel to a 
neighbouring voxel at a mean rate of d where d = ^. In other words, within an infinitesimal time 
At, the probability that a signalling molecule diffuses to a neighbouring voxel is d At. Note that 
the expression for d can be derived from spatially discretising the diffusion equation into regular 
cubic voxels of volume W'^, see [HI p.341] or [IHl Section 3]. 

The dashed lines in Figure [^indicate the boundary of our transmission medium. Many different 
boundary conditions are used by engineers and physicists to model what happens when a molecule 
reaches the boundary of a medium. Two typical boundary conditions are absorbing and reflecting 
boundaries BH. An absorbing boundary means that a molecule can leave the transmission medium 
and once the molecule has left, it will not return to the medium. For example, in Figure we 
allow molecules to leave the medium via one surface of Voxel 16 as indicated by the one-way arrow. 
Mathematically, this is modelled by a rate of leaving the medium, similar to that of modelling 
the diffusion between the voxels. A reflecting boundary means that a molecule cannot leave the 
medium, i.e. a molecule hitting a reflecting boundary will stay in the voxel. Our model can capture 
these boundary conditions. 

It has been shown in [38l [31] that in order for RDME to produce physically meaningful results, 
the voxel dimension W cannot be too small and in order for RDME to reduce the discretisation 
error, W cannot be too large. In this paper, we assume that W comes from a valid range. The choice 
of W is beyond the scope of the paper and the reader can refer to |3Sl El] for further discussion. 

For simplicity, we assume that the medium is homogeneous with a constant diffusion coefficient 
D. It is straightforward to extend the framework to cover inhomogeneous medium [35]. It is also 
possible to use non-cubic voxels, see [HI IT3] . 

3.1.2 Transmitter 

We assume the transmitter occupies one voxel. However, it is straightforward to generalise to the 
case where a transmitter occupies multiple voxels. We limit our consideration to one symbol interval 
without inter-symbol interference (ISI) at this moment. We will discuss the multiple symbol interval 
case with ISI in Section llTSl 

We assume that the transmitter can send K different symbols s = 0, l,..,iF — 1 where each 
symbol s is characterised by an emission pattern Us(t). The role of emission pattern in molecular 
communication is the same as that of transmitted signal in electromagnetic communication. If a 
transmitter uses a deterministic emission pattern Us{t) to represent symbol s, it means the trans- 




mitter emits Us{t) signalling molecnles into the transmitter voxel at time t. We nse an example to 
illnstrate the meaning of emission pattern. Consider an emission pattern ui(t) for Symbol 1 where 
ui(t) = 6t^i.2 + (5i,5.6 + 2(5t,8.i where denotes the Kronecker delta, which means 6ij = 1 if and only 
ii i = j. The emission pattern Ui{t) means that, for Symbol 1, the transmitter emits one signalling 
molecnle at times 1.2 and 5.6, two signalling molecnles at time 8.1 and does not emit any molecnles 
at any other times. 

In this paper, we assnme that the emission pattern for each symbol is prodnced by a set of 
chemical reactions located in the transmitter voxel. Given K symbols, the transmitter nses K 
different reactions to generate these symbols, see Figure As an example, a class of chemical 
reactions inside living cells |16] is 

RNA ^ RNA + A (1) 

where ribonucleic acid (RNA, which is a molecule commonly found in living cells) produces the 
chemical species A. This class of chemical reactions can be modelled by a Poisson process where 
molecules of A are produced at a mean rate of k m- Note that the emission patterns produced 
by chemical reactions are not deterministic, but stochastic. The mean emission pattern of this 
chemical reaction is E[n(t)] = k. 

Following on from the above example, one can realise Amplitude Shift Keying (ASK) in molec¬ 
ular communication by using different chemical reactions that can produce signalling molecules at 
different mean rates. For example, if there are four different reactions that can produce signalling 
molecules at four different mean rates of ko, ki, K 2 and ks, then one can use these four different 
reactions to produce 4 different symbols. Note that it is possible for the four chemical reactions to 
produce the same emission pattern (or realisation), though with different probabilities. 

A standard result in physical chemistry shows that the dynamics of a set of chemical reactions 
can be modelled by a CTMP [18]. Therefore, we will model the transmitter by a CTMP. Note 
that, in this paper, we will not specify the sets of chemical reactions used by the transmitter except 
for simulation because the MAP demodulator does not explicitly depend on the sets of chemical 
reactions that the transmitter uses. 

3.1.3 Receiver 

We assume the receiver occupies one voxel and we use R to denote the index of the voxel at which 
the receiver is located. In Figure we assume the receiver is at Voxel 7 (light grey) and hence 
R = 7 for this example. In addition, we assume that the transmitter and receiver voxels are distinct. 
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We assume that the receiver has M non-interacting receptors and we use E as the chemical 
name for an unbound receptor. These receptors are hxed in space and do not diffuse, and they 
are only found in the receiver voxel. Furthermore, these receptors are assumed to be uniformly 
distributed in the receiver voxel. 

The receptor E can bind to a signalling molecule S to form a ligand-receptor complex (or complex 
for short) C, which is a molecule formed by combining E and S. This is known as ligand-receptor 
binding in molecular biology literature |19]. The binding reaction can be written as the chemical 
equation: 

S + E 4 C (2) 


where A is the reaction rate constant. Since the receptors are only found in the receiver voxel, the 
binding reaction occurs in a volume of W^, which is the volume of a voxel. The rate at which the 
complexes C is formed is given by the product of the number of signalling molecules in the 
receiver voxel and the number of unbound receptor^ We dehne A = and will use A in the 


CTMP. Note that this is equivalent to ligand-receptor binding model used in [301 Section V-B]. 

A ligand-receptor complex C can dissociate into an unbound receptor E and a signalling molecule 
S. This can be represented by the chemical equation 


C AE-fS 


(3) 


where /i is the reaction rate constant. The rate at which the complexes are dissociating is given by 
the product of /i and the number of complexe^ 

Since a receptor can either be in an unbound state E or in a complex C, we have the following 
conservation relation: the number of unbound receptors plus the number of complexes is equal to 
the total number of receptors M. 

^ This footnote explains how comes about. Consider a chemical reaction where reactants S and E react to 
form product C. We assume the reactions are taking place within a volume of W^. Let cs, ce and cc be, respectively, 
the concentration of S, E and C in the volume. The law of mass action says that = A ce cs. In the case of the 
CTMP or RDME in this paper, we want to keep track of the number of molecules in a volume (the voxel) instead. 
Let ns, riE and nc be, respectively, the number of S, E and C molecules in the volume. Since concentration and 
molecule counts are related by ccW^ = nc etc, we will, in a mathematically loose way, write ns n^. 

Since uq is a discrete quantity, the derivative is not defined but we can interpret it as the production rate of 
molecules C. This explains how to convert the law of mass action, which is in terms of concentration, to the rate 
law used in RDME which is in terms molecular counts. This conversion is also discussed in [sniiis]- 

^The law of mass action for the dissociation reaction is = —/i cc where cc is the concentration of the 
complexes. We can use the same argument in Eootnotej^to show that the dissociation rate of C is /i nc where nc 
is the number of complexes. In particular, note that no scaling by volume of the voxel is required. 
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3.2 General end-to-end model 


In order to derive the MAP demodulator, we need an end-to-end model which includes the trans¬ 
mitter, the medium and the ligand-binding process, see Figure Since chemical reactions (which 
includes the chemical reactions in the transmitter as well as the ligand-receptor binding process in 
the receiver) and diffusion can be modelled by CTMP, it is possible to use a CTMP as an end-to- 
end model. In this section we present a general end-to-end model that includes the transmitter, 
diffusion and the ligand-receptor process in the receiver. An excellent tutorial introduction to the 
modelling of chemical reactions and diffusion by using CTMP can be found in [13]. We have also 
included an example in Appendix 

The aim of the end-to-end model is to determine the properties of the receiver signal from the 
transmitter signal. The receiver signal in our case is the number of complexes over time. Since the 
transmitter uses K symbols, the transmitter signal is generated by one of the K sets of chemical 
reactions. This means that we need K end-to-end models with a model for each of the K symbols or 
sets of chemical reactions. The principle behind building these K models is identical so without loss 
of generality, we will assume that the model here is for Symbol 0. We begin with a few definitions. 

Let n.i{t) (where 1 ^ i ^ Ny) he the number of signalling molecules S in Voxel i at time t. In 
particular, since we have defined R to be the index of the receiver voxel, rinit) is the number of 
signalling molecules in the receiver voxel. We assume the transmitter is a set of chemical reactions 
which uses H intermediate chemical species Qi, Q 2 , ... and Qh and these intermediate species 
remain in the transmitter voxel. Let nQ.{t) be the number of chemical species Qi in the transmitter 
voxel at time t. Molecules may also be degraded or leave the system forever if absorbing boundary 
condition is used. We use n^(t) to denote the cumulative number of molecules that have left the 
system. Note that since ni(t),nQ.(t) and n^(t) are molecular counts, they must belong to the set 
of non-negative integers Z^o- We define the vector N{t) e to be: 


N{t) 


ni{t) ... n 7 v„(t) nQ^{t) 


nQnit) 


UAit) 


T 


(4) 


where ^ denotes matrix transpose. 

Let b{t) denote the number of complexes or bound receptors at time t and Z^^m] denote the set 
of integers between 0 and M inclusively. We require that a valid b{t) must be an element of Zp^^]. 
Since there are M receptors, the number of unbound receptor is M — b{t). 

The state of the end-to-end model is the tuple {N{t),b(t)). We will now specify the transition 
probabilities from state {N(t), b{t)) to state {N(t + At), b{t+At)). State transitions can be caused by 
any one of these events: a chemical reaction in the transmitter, the diffusion of a signalling molecule 
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from a voxel to neighbouring voxel, and the binding or unbinding of a receptor in the receiver. We 
know from the theory of CTMP that the probability of two events taking place in an inhnitesimal 
duration of At is of the order of o(At^). Intuitively, this means only one event can occur within At. 
We can divide the transition probabilities from {N(t), b{t)) to {N(t + At), b{t + At)) into 2 groups 
depending on whether the number of complexes has changed or not in the time interval (t, t + At). 
If the number of complexes has changed from time t to t + At, i.e. b(t + At) ^ b(t), this means 
either a binding reaction (|^ or a unbinding reaction (|^ has occurred. 

If a binding reaction (|^ has occurred, then the number of signalling molecules in the receiver 
voxel is decreased by 1 and the number of complexes b(t) is increased by 1. This reaction occurs at 
a mean rate of A nuit) (M — b(t)). We use 1* to denote the standard basis vector with a ‘1’ at the 
i-th position. We can write the state transition probability of the receptor binding reaction ([^ as: 

P[N{t + At) = N{t) - 1 r, bit + At) = b{t) + l|iV(t), b{t)] = A nR{t) (M - b{t)) At (5) 

Recalling that R is the index of the receiver voxel and njiit) is the R-th element of iV(t) in Q, 
the expression A(t + At) = Nit) — 1r is equivalent to riRit + At) = nuit) — 1, which means the 
number of signalling molecules in the receiver voxel has decreased by 1. Similarly, the expression 
6(t + At) = 6(t) + 1 says the number of complexes has increased by 1. The right-hand side (RHS) 
of (|^ is the transition probability and is given by the product of mean reaction rate and At. 

Similarly, the transition probability of the unbinding reaction is given by: 

P[A(t -I- At) = Nit) + Ir, bit + At) = bit) — 1| A(t), bit)] = p 6(t) At (6) 

where RHS of (|^ is the transition probability. 

We now specify the second group of transition probabilities with 6(t -I- At) = 6(t). These 
transitions are caused by either a reaction in the transmitter or diffusion of signalling molecules 
between neighbouring voxels. Let rii,rij e be two valid Nit) vectors; let also (3 e 

For rji ^ rjj, we write 

P[A(t -f At) = rji, bit + At) = /3\Nit) = rjj, bit) = /3] = dij At (7) 

where dij is the transition rate from state (?7j,/d) to state (r;,,/?). Since this transition is due to 
either a reaction in the transmitter or diffusion, dij is independent of the number of complexes (3. 
Depending on the type of transition, the value of dij can depend on the reaction constants in the 
transmitter, diffusion rate and some states of rjj. For example, if the transition from r]j to rji is 
caused by the diffusion of a signalling molecule from Voxel 1 to Voxel 2, we have rji = rjj — li +12 at 
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a rate of drjj^i where rjj^i is the first element in rjj or equivalently the number of signalling molecules 
in Voxel 1 in state rjj] so, for this example, dij = drjj^i. The main advantage of using Equation ([^ 
is that it allows us a cleaner abstraction to solve the Bayesian hltering problem when deriving the 
MAP demodulator. We also remark that we will not specify the exact expression of dij because 
dij^s do not appear explicitly in the demodulator. 

Equations ([^, ([^ and Q specify all the possible state transitions. The probability of no state 
transition is therefore: 


P[iV(f + At) = r]j,b(t + At) = b{t)\N(t) = rij,b{t)] 

= 1 — djj At — X nnit) (M — b{t)) At — ^ b{t) At (8) 

where 

(^) 

i^j 

We have now specified all the state transition probabilities for Symbol 0. If a different symbol 
is used, the value of H, the dimension of N{t) and the dij parameters can change. However, the 
state transition probabilities still can be summarised by Equations of the form ([^, ([^ and ([^. In 
any case, the derivation of the MAP demodulator only requires us to work with one symbol at a 
time. Hence, we will use Equations (|^-([^ for any transmission symbol. 

3.3 Discussion 

Note that the CTMP includes all the three sources of noise in our system, due to chemical reactions 
in the transmitter, random diffusive movements in the medium and the ligand-receptor binding 
process at the receiver. Some of these noise components have also been characterised in earlier 
literature. For example, [52] discusses sampling noise at the transmitter and counting noise at the 
receiver. One can study these noises using the derived CTMP and let us take sampling noise as an 
example. For the moment, let us isolate the transmitter voxel from the propagation medium, i.e. 
we do not allow the signalling molecules to leave or enter the transmitter voxel. In this case, the 
number of signalling molecules in the transmitter voxel is due entirely to chemical reactions and 
we use nigoiated(^) to denote the number of signalling molecules in the isolated transmitter voxel at 
time t. Now, let us consider the transmitter voxel again but we allow signalling molecules to diffuse 
in and out of the voxel; we use Uconnected(t) to denote the number of signalling molecules in the 
transmitter voxel in this case. It is natural to consider Uconnected(^) as the transmitter signal because 
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this is the number of signalling molecules in the transmitter voxel. However, /iconnected(^) can be 
different from nisoiated(^) because signalling molecules can diffuse in and out of the transmitter voxel, 
which is the cause of sampling noise. 

Equations (|^ to Q hold for any valid state b{t)). If we collect all the transition probability 

equations for all valid states, then we can form the inhnitesimal generator of the CTMP. For a given 
initial probability distribution of the initial state (A^(0), 6(0)), one can in principle solve the hrst 
order ordinary differential equation (ODE) associated with the inhnitesimal generator to compute 
the probability of the number of complexes b(t), or the property of the receiver signal. However, 
in practice, this ODE is of a very high dimension and it is an active area of research to derive 
algorithms to solve this ODE efficiently and accurately [53] . 

4 The MAP demodulator 

This section aims to derive the optimal demodulator using the CTMP derived in the previous 
section. We assume the input to the demodulator is the continuous-time signal b{t) which is the 
number of complexes at time t. There are a number of reasons why we choose to work with the 
continuous-time signal b{t), rather than its sampled version. First, the signal b{t) may not be strictly 
band limited in the frequency domain. Second, our results show that the optimal demodulator needs 
to know the time instances at which a receptor is switching from the unbound to bound state. This 
timing information, which is essentially an impulse, is unfortunately lost by sampling b{t). Third, 
the solution of the proposed decoding problem can be used to benchmark molecular circuit [29] based 
decoders. Since molecular circuits use chemical reactions for computation, they are fundamentally 
analogue circuits. Fourth, there is an increasing interest in the circuit design community to design 
low-power analogue signal processing circuits iai. 

An intermediate step to derive the MAP demodulator is to solve a continuous-time hltering 
problem. In a hltering problem, one uses all the observations available up till time t to predict 
the system state at time t. For our case, the observations are the number of complexes b{t) in 
the continuous interval [0,t]. (This means the number of observations is inhnite because we are 
considering the continuous-time signal b(t) in a non-zero time interval [0,t].) We use B{t) = 
{6(r); 0 ^ T ^ t} to denote the section of the signal b(t) in the time interval [0, t]. (Note that the {} 
is not a set notation. Here B{t) is a realisation of the number of complexes in [0,f] of the CTMP. 
This notation is used in some non-linear hltering literature, e.g. [55].) We can think of B{t) as the 
history of the number of complexes up till time t. The demodulation problem is to use the history 
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B{t) to determine which symbol the transmitter has sent. 


In Sections |4.1| to |4dll we consider only one symbol interval and do not consider ISI. We consider 
the ISI case in Section llThl 


4.1 The MAP framework 

We adopt a MAP framework for detection. Let P[s|i3(t)] denote the posteriori probability that 
symbol s has been sent given the history B{t). If the demodulation decision is to be done at time 
t, then the demodulator decides that symbol s has been sent if 


s = argmax,^o_ 

Instead of working with P[s|I3(t)], we will work with its logarithm. Let 

L,{t) = log(P[s|S(t)]) 


( 10 ) 


( 11 ) 


The first step is to determine Ls{t + At) from Lgit). Given B{t) is the section of h{t) in the time 
interval [0, t], one can consider B{t + At) as the concatenation of B{t) and the section of b(t) in the 
time interval (t, t + At]. Over an inhnitesimal At, we can consider the signal h{t) is a constant in 
the time interval (t, t + At]; we therefore abuse the notation and use h{t + At) to denote the section 
of 6(t) in the time interval (t, t + At]. By using Bayes’ rule, it can be shown that 

Ls{t + At) =Ls{t) + log(P[6(t + At)|s, B{t)]) — log(P[6(t + At)|i3(t)]) (12) 


where P[&(t + At)\s, B{t)] is the probability that there are b{t + At) complexes given that the 
transmitter has sent the symbol s and the previous history B{t). The last term on the RHS of 


(12), i.e. P[6(t + At)|i3(t)], is independent of the transmission symbol so we do not need it for the 
purpose of detection. We will focus on determining P[6(t + At)|s, B{t)). 


4.2 Computing P[b{t + At)|s, B{t)] 

The problem of determining the probability P[&(t + At)\s, B{t)] is essentially a Bayesian hltering 
problem. Recall that the complete state of the system is {N(t),b(t)) and the receiver can only 
observe b(t), therefore the task of the receiver is to use the history B{t) and the system model to do 
prediction. Standard method can be used to derive P[6(t + At)\s,B{t)] but the derivation is long. 
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especially because of the diffusion terms; the derivation can be found in Appendix The result is 


P[b{t + At)\s,B{t)] 

=<^6(t+Ap,b(i)+i A(M — b{t)) At Ei[nR{t)\s, B{t)] + 5b{t+At),b{t)-i f^b{t) At + 
Sb{t+At)Mt) (1 “ “ b{t))E[nR{t)\s, B{t)] At - fib{t) At) 


(13) 


Note that only one of the three terms on the RHS of Equation ( |13[ ) is non-zero depending on whether 
the observed b{t + At) is one more, one less or equal to that of b(t). The term E[nj:j(t)|s, B{t)] is the 
expected number of signalling molecules in the receiver voxel given the history and the symbol s. 
The meaning of this term is that the receiver uses the history to predict what the expected number 
of signalling molecules in the receiver voxel is. Note that only the chemical kinetic parameters A and 


fi of the receptor appear explicitly in Equation (13). Other parameters, such as the set of chemical 


reaction that generate Symbol s and the diffusion coefficient, do not appear explicitly in Equation 


(13) but influence the system behaviour via the term E[nij(t)|s, i3(t)]. 


4.3 The demodulation filter 


By substituting Equation (13) into Equation (12) and let At go to zero, we show in Appendix [b| 
that 


dL,{t) dU{t) 


log(E[nK(t)|s,i3(t)]) - A(M - b{t))E[nR{t)\s,B{t)] + L{t) 


(14) 


dt dt 

with Ls(0) initialised to the logarithm of the prior probability that Symbol s is sent. Equation 


(14) is the optimal demodulation hlter. The term U{t) is the cumulative number of times that 


the receptors have turned from the unbound to bound state. The meaning of U{t) is illustrated 
in Figure assuming there are two receptors. The top two pictures in Figure show the state 
transitions for the two receptors. The third picture shows the function U(t) which is increased by 
one every time a receptor switches from the unbound to bound state. The bottom picture shows 
which is the derivative of U{t). Note that consists of a train of impulses (or Dirac deltas) 
where the timings of the impulses are the times at which a receptor binding event occurs. Loosely 
speaking, one may also view as max(^^^, 0). 


The function L{t), which is the last term on the RHS of (14), contains all the terms that are 


independent of Symbol s. Since Ls{t) does not appear on the RHS of (14), this means that L{t) 
adds the same contribution to all Ls(t) for all s = 0,K — 1. We can therefore ignore L(t) for the 
purpose of demodulation. 
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The term E[n/^(t)|s, i3(t)] in Equation (14) is the prediction of the mean number of signalling 
molecules in the receiver voxel using the history of receptor state. This is a hltering problem which 
requires extensive computation. Instead, we assume that the receiver has prior knowledge that 
if Symbol s is transmitted, then the mean number of signalling molecules in the receiver voxel is 
o's{t) = E[?7,j^(t)|s] and the receiver uses (Js(t) for demodulation. We can view Csit) as internal 
models that the demodulator uses. The use of internal models is fairly common in signal processing 
and communication, e.g. a matched filter correlates the measured data with an expected response. 

After making the modifications described in the last two paragraphs, we are now ready to 
describe the demodulator. Using b(t) as the input, the demodulator runs the following K continuous¬ 
time hlters in parallel: 

dZs{t) dU{t) 


dt 


dt 


log(tr,(«)) - \(M - b{t))a,(t) 


( 16 ) 


where ^^(0) is initialised to the logarithm of the prior probability that the transmitter sends Symbol 
s. If the demodulator makes the decision at time t, then the demodulator decides that Symbol s 
has been transmitted if 


s = argmax^^o 


(16) 


The demodulator structure is illustrated in Figure]^ By comparing Equations (14) and (15), it can 
be shown that Lg^^it) — Ls 2 it) = Zg^it) — Zg 2 {t) for any two symbols si and S 2 - An interpretation of 
the demodulation filter output Zg(t) is that exp(Zs(f)) is proportional to the posteriori probability 
P[.|i3(t)]. 


Note that the replacement of E[?7,K(f)|s, S(f)] in Equation (14) by ag{t) means that the de¬ 


modulation filter (15) is sub-optimal. If E[nH(t)|s,-B(f)] and ag(t) are close to each other, we 


expect the performance degradation to be small. It is an open problem what the difference be¬ 
tween E[?7,7j(f)|s, B{t)\ and ag{t) is. The difficulty in answering this problem is due to the fact that 
ligand-receptor binding is a nonlinear process. In Appendix we motivate the closeness between 
E[nij(t)|s, B{t)\ and ag{t) by studying an analogous problem in linear time-invariant (LTI) systems. 
We will compare the performance of the optimal and sub-optimal demodulation filters in Section 

o 


In order to understand Equation (15), we consider the situation where Symbol 1 generates a lot 


more signalling molecules than Symbol 0 such that it results in more signalling molecules in the 
receiver voxel, or (Ti{t) > ao(t) for all t. If the transmitter sends Symbol 1, then more signalling 
molecules are expected to reach the receiver voxel. The consequence is that there are more receptor 
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binding events and the number of unbound receptors (M — b{t)) is smaller. Therefore, in Equation 


(15), we expect a big positive contribution from the hrst term on the RHS and a small negative 
contribution from the second term. The net effect is a big Zi{t). On the other hand, if the 
transmitter sends Symbol 0, the number of receptor binding events is smaller and (M — h{t)) is 
big. This results in a smaller ZQ(t). Therefore, Zi{t) is likely to be bigger than Z^it), which means 
correct detection. 


4.4 Discussions 


4.4.1 Implementation issues 

The implementation of the demodulator is an open research problem. The demodulation hlter 


(15) is an analogue hlter and it requires the internal model crj(t). The design of analogue circuits 

for a recent overview. 


using chemical reactions for calculations is an active research area, see 
The demodulation hlter requires logarithm, multiplication, subtraction, integration and counting 
the number of times the receptors have switched from the unbound to bound state. An analogue 
molecular circuit for calculating logarithm is presented in 0. There are also circuits that can 
perform multiplication and subtraction [SS]- Living cells are known to use integration [SB]. It 
may be possible to implement the counting using a chemical reaction [39] . It may be possible to 
approximate the internal models by using some lower order chemical reactions. This discussion 
shows that some components to implement the demodulator exist but the exact implementation 
remains an open problem. 

In order to bypass the difficulty of computing E[nij(t)|s, we have proposed to use internal 

models <Js(t). An open research problem is to study sub-optimal estimation of E[nij(f)|s, 

Note that |BS| shows that analogue computation is more efficient than digital computation if high 
precision is not required. An interesting problem is to study the impact of low precision analogue 
calculations on the demodulation performance. 


4.4.2 Continuous transmission medium versus voxels 

We mention in Section [^ that some papers in molecular communications assume a continuous 
medium (rather than discretising the medium into voxels) and a non-zero receiver size. If a con¬ 
tinuous medium is assumed, the state of a signalling molecule in the transmission medium is its 
position. Let p{x, t) be the probability that a molecule is at position x at time t. The time evolution 
of p{x, t) can again be modelled by a CTMP, which is continuous in both time and space. The time 
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evolution of p{x,t) is described by the differential Chapman-Kolmogorov Equation (CKE) [H]. It 
is possible to derive an alternative CTMP by replacing the diffusion of signalling molecules in Equa¬ 
tion Q by differential CKE as well as by adding equations to describe how the signalling molecules 
enter or leave the receiver voxel. We can then apply Bayesian hltering to this alternative CTMP. We 
expect this alternative CTMP will give the same demodulator hlter because in our derivation based 
on discrete voxels, the MAP demodulator depends only on the number of the signalling molecules in 
the receiver voxel and the diffusion parameters do not appear explicitly in the demodulation hlter. 
Further support of this argument is given in the derivation in Appendix where we show that the 
dij parameters in Q are cancelled out in deriving the Bayesian hlter. 


4.5 Inter-symbol interference 


It is in principle possible to use the demodulation hiters (15) to deal with the case with ISI. We 


use Tj, to denote one symbol duration and we assume that ISI only lasts for a hnite amount of 
time. Specihcally, consider a symbol sent in [kTx,{k + l)Tx], we assume that the ehect of this 
transmission can be neglected after the time {k + n)Tx. This can be realised by appropriately 
choosing the transmitter and receiver parameters, and n. 

In order to make the explanation here a bit more concrete, we assume that the transmitter 
uses K = 2 symbols and n = 3. Over a duration of n = 3 symbols, the possible sequences 
sent by the transmitter are 000, 001, 010, 011, ... , 111. Let (Jofifiit) denote the mean number 
of signalling molecules at the receiver voxel if the sequence 000 is sent. We can similarly dehne 
o'o,o,i(^), ••• Consider the transmission of three consecutive symbols Sk- 2 ,Sk-i and Sk- As¬ 

suming that we have an estimation of the hrst two symbols Sk -2 and Sk-i, then the decoding of Sk 


can be done by using the demodulation hlter (15) by replacing as{t) by o's^._ 2 ,sk-i,s- For example, 
if Sk -2 = 1 and = 0, then one can decode what Sk is by using the demodulator hiters cTi_o,i(t) 
and Although the decision feedback based method can solve the ISI problem, the number 

of internal models increases exponentially with the memory length parameter n. 

The reason why we need to consider all 2"’ possible transmission sequences is that the ligand- 
receptor binding process has a non-linear reaction rate. A method to reduce the number of internal 
models is to design the system so that cto,o,o(^) etc. can be decomposed into a sum. Let as(t) 
(s = 0,1) be the mean number of signalling molecules at the receiver voxel if the symbol s is sent 
for one symbol duration and in the absence of ISI. If 


*^S1,S2,S3 (f) ^ -1- (Jg^it Tx) T ^33 (^) 


(17) 
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holds for all Si,S 2 and S 3 , then one can again make nse of decision feedback to decode the ISI 


signal. However, this time, only K internal models are needed. Eqnation (17) can be made to 


hold approximately if the nnmber of receptors is large. This can be explained as follows. First 
of all, if ligand-receptor binding is absent, this means there is only free diffnsion then Eqnations 


(17) holds becanse the mean nnmber of signalling molecnles obeys the diffnsion eqnation which is 


linear. This means that we need to create an environment that “looks like” free diffnsion even when 
ligand-receptor binding is present. This can be realised if the nnmber of signalling molecnles that 
are bonnd to the receptors is small compared to those that are free. A method to achieve this is to 
increase the nnmber of receptors. We will demonstrate this with a numerical example in Section 
However, it is still an open problem to solve the ISI in the general case. 

5 Properties of the demodulator 

The aim of this section is to study the properties of the MAP demodulator numerically. We begin 
with the methodology. 


5.1 Methodology 


We assume the diffusion coefficient D of the medium is 1 /im^s“^. The receptor parameters are A 
= 0.005 s“^, A = and fi = I s“^. These values are similar to those used in [50] and [30] 

The above parameter values will be used for all the numerical experiments. 

For each experiment, the transmitter uses either K = 2 or K = 3 symbols. Each symbol is 
generated by a different sets of chemical reactions. Different experiments may use different sets of 
chemical reactions and will be described later. The number of receptors also varies between the 
experiments. 

We use the Stochastic Simulation Algorithm (SSA) [59] to obtain realisations of b{t) which is 
the number of complexes over time. SSA is a standard algorithm in chemistry to simulate diffusion 
and reactions; it is essentially an algorithm to simulate a CTMP. 


In order to use Equation (15), we require the mean number of signalling molecules as{t) in the 


receiver voxel when Symbol s is sent. Unfortunately, it is not possible to analytically compute 


^ The paper [80] considers ligand-receptor binding in the chemical master equation setting. In our notation, the 
parameter values in [3^ are D = 100 /im^s“^, A = 0.2 /rm^ s“^ and /i = 10 s“^. These parameters are 10-100 times 


faster than ours and can be considered as a time-scaling. Note that 
and /i. 


uses k+ and instead of, respectively, A 
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(Js(t) from the CTMP because of moment closure problem which arises when the transition rate is 
a non-linear function of the state |6^. We therefore resort to simulation to estimate Csit). Each 
time when we need an (Ts(t), we run SSA simulation 500 times and average the results to obtain 
(Ts{t). Note that these simulations are different from those that we use to generate b(t) for the 
performance study. In other words, the simulations for estimating as{t) and for performance study 
are completely independent. 

Once b{t) and as(t) are obtained, we use numerical integration to calculate Zs(t) using Equation 


(15). We assume that all symbols appear with equal probability, so we initialise Zg{Q) = 0 for all s. 


5.2 Optimal filter (14) versus sub-optimal filter (15) 


The optimal demodulation hlter (14) requires the term E[nj?(t)|s, i3(t)] which can be obtained by 


solving a Bayesian hltering problem. Since hltering problems are computationally expensive to 


solve, we propose the sub-optimal demodulation hlter (15) which uses as(t) = E[?7,jj(f)|s] as an 


internal model. The aim of this section is to compare the performance of these two demodulation 
hlters. 

In this comparison, we consider a medium of l/rm x l/rm x l/rm. We assume a voxel size of 
(|/im)^ (i.e. IT = I /^m), creating an array of 3 x 1 x 1 voxels. The voxel co-ordinates of transmitter 
and receiver are, respectively, (1,1,1) and (3,1,1). A rehecting boundary condition is assumed. 

The reason why we have chosen to use such a small number of voxels is because of the dimen¬ 
sionality of the hltering problem. For example, if each voxel can have a maximum of 100 signalling 
molecules at a time, then there are 10® possible N{t) vectors and the hltering problem has to 
estimate the probability P[A^(t)|s, i3(f)] for each possible N{t) vector. Although there are approx¬ 
imation techniques to solve the Bayesian hltering problem, that would introduce inaccuracies. The 
use of small number of voxels will allow us to compute P[n/j(f)|s, B{t)\ precisely. 

For this experiment, we use K = 2 symbols and two values of M (the number of receptors): 5 
and 10. Both Symbols 0 and 1 use Reaction ([^ such that Symbols 0 and 1 causes, respectively, 
10 and 50 signalling molecules to be generated per second on average by the transmitter. The 
simulation time is about 1.8 seconds. 

We hrst show that E[nij(f)|s, i3(f)] and as{t) are rather similar. Figure shows as{t) and a 
trajectory of E[n/{(f)|s, 13(f)] (obtained from one realisation of b(t)) are pretty similar. This result 
is obtained from using M = 10 and Symbol 1. The results for other choices of M, transmission 
symbols or other realisations of b{t) are similar. 
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Figure 1^ shows the mean symbol error rates (SERs), for both optimal and sub-optimal demod¬ 
ulation hlters, if the detection is done at time t = 1, 1.05, 1.1, ..., 1.8. The SERs is obtained from 
400 realisations of b(t). The difference in SERs between the optimal and sub-optimal hlter is less 
than 1%. We have also checked that the two demodulators make the same decoding decision on 
average 99.3% of the time. 


In the rest of this section, we will use the sub-optimal demodulation hlter (15) because of its 
lower computational complexity. 


5.3 Properties of the demodulator output 

We consider a medium of 2/rm x 2/im x 1 /rm. We assume a voxel size of (|pm)^ (i.e. IF = | /^rn), 
creating an array of 6 x 6 x 3 voxels. The transmitter and receiver are located at (0.5,0.8,0.5) and 
(1.5,0.8,0.5) (in /rm) in the medium. The voxel co-ordinates are (2,3,2) and (5,3,2) respectively. We 
assume an absorbing boundary for the medium and the signalling molecules escape from a boundary 
voxel surface at a rate of This conhguration will be used for the rest of this section. 

For this experiment, we use K = 2 symbols and M = 50 receptors. Both Symbols 0 and 1 use 
Reaction ([^ such that Symbols 0 and 1 causes, respectively, 40 and 80 signalling molecules to be 
generated per second on average by the transmitter. The simulation time is about 3 seconds. 

Figure [^shows the demodulation hlter outputs Zo{t) and Zi{t) if the transmitter sends a Symbol 
0. It can be seen that Zo{t) > Zi(t) most of the time after t = 1.2, which means the detection is 
likely to be correct after this time. The sawtooth like appearance of Zo{t) and Zi(t) is due to the 
fact that every time when a receptor is bound, there is a jump in the hlter output according to 
Equation ( [l5| . Figure [^shows the hlter outputs Zo{t) and Zi{t) if the transmitter sends a Symbol 
1; the behaviour is similar. 

Figure l^shows the mean hlter outputs ZQ{t) and Zi{t) if the transmitter sends a Symbol 0. The 
mean is computed over 200 realisations of b{t). It can be seen that the mean hlter output of Zo(t) 
is greater than that of Zi{t). Similarly, if Symbol 1 is sent, then we expect of the mean of Zi{t) to 
be bigger. The hgure is not shown for brevity. 

Figure [TO] shows the mean SERs for Symbols 0 and 1 if the detection is done at time t. The SER 
for Symbol 1 is high initially but as more information is processed over time, the SER drops to a 


low value. This experiment shows that it is possible to use the analogue demodulation hlter (15) 


to compute a quantity that allows us to distinguish between two emission patterns at the receiver. 
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5.4 Impact of number of receptors 


We continue with the setting of |5.3| but we vary the number of receptors between 1 and 20. We 
assume the demodulator makes the decision ait = 2.5 and calculate the mean SER for both symbols 
at t = 2.5. Figure plots the SERs versus the number of receptors. It can be seen that the SER 
drops with increasing number of receptors. 

We have used K = 2 symbols so far. We retain the current Symbols 0 and 1, and add a Symbol 
2 which is also of the form of Reaction ([^ but its mean rate of production of signalling molecules 
is 3 times that of Symbol 0. The number of receptors M used are: 1, 10, 20, ..., 150. We compute 
the average SER at t = 2.5 assuming each symbol is transmitted with equal probability. We plot 


the logarithm of the average SER against log(M) in Figure 12 It can be seen that the SER drops 
with increasing number of receptors M. The plot in Figure [T^ suggests that, when the number of 
receptors M is large, the relationship between logarithm of SER and log(M) is linear. We perform 
a least-squares fit for M between 50 and 150. The fitted straight line is shown in Figure [T^ and it 
has a slope of —1.13. A possible explanation is that, because the receptors are non-interacting, each 
receptor provides an independent observation. The empirical evidence suggests that the average 
SER scales according to A asymptotically provided that the voxel volume can contain that many 
receptors. 


5.5 Distinguishability of different chemical reactions 

Equation suggests that if the transmitter uses two sets of reactions which have almost the same 
mean number of signalling molecules in the receiver voxel, then it may be difficult to distinguish 
between these two symbols. In this study. Symbol 0 is generated by Reaction Q with a rate of k 
while Symbol 1 is generated by: 


[RNA]o„ « [RNAJopp 

(18) 

[RNA]o^ i [RNAJo,, + S 

(19) 


where we assume that RNA can be in an ON or OFF state, and signalling molecules S are only 
produced when the RNA is in the ON-state. We assume that the there is an equal probability for the 
RNA to be in the two states and the reaction rate constant for the production of signalling molecule 
S from [RNAJqj.^ is 2 k. This means that the mean rate of production of signalling molecules S by 


Symbols 0 and 1 are the same. This gives rise to very similar cro(f) and cTi(f). Figure 13 shows the 
demodulation filter outputs Zo{t) and Zi{t) for one simulation. It can be seen that the two outputs 
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are almost indistinguishable. Consequently, the SER is pretty high. This shows that symbols 
generated by reactions which have similar mean number of signalling molecules at the receiver 
voxel can be hard to distinguish. 


5.6 Inter-symbol interference 

The aim of this experiment is to study the performance of the demodulator in the presence of ISI. 
We use the decision feedback method described in Section |4.5| together with the approximation 
decomposition in ([T^. We vary the number of receptors from 25 to 150. We use two different 
memory lengths £ of 4 and 5. If the memory length is £, we express the mean number of output 


molecules at the current symbol interval as a sum of £ terms, i.e. a generalisation of (17) to £ terms. 
We carry out the simulation for a duration of 20 symbol lengths and compute the average SER over 


20 symbols. Figure 14 shows the average SER versus the number of receptors. It can be seen that 
an increasing number of receptors can also be used to deal with ISI. 


6 Conclusions and future work 

This paper studies a diffusion-based molecular communication network that uses different sets of 
chemical reactions to represent different transmission symbols. We focus on the demodulation 
problem. We assume the receiver uses a ligand-receptor binding process and uses the continuous 
history of the number of ligand-receptor complexes over time as the input signal to the demodulator. 
We derive the maximum a posteriori demodulator by solving a Bayesian hltering problem. 


A Proof of Equation ( 13 ) 


Let s denote the transmitted symbol, our aim is to determine P[6(f -I- At)\s,B{t)] in terms of the 
quantity at time t. Recalling that {N{t),b(t)) is the state of the CTMP and since only B{t) is 
observed, the problem of predicting h{t + At) from B{t) is a Bayesian hltering or hidden Markov 
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model problem. The first step is to condition on the state of the system, as follows: 

P[b{t + At)\s,B{t)] (20) 

= J]P[N{t + At)=r],,bit + At)\s,B{t)] (21) 

i 

= Y,Y,^[N{t + At)=1^,Mt + ^t)\s,N{t)=1^j,B{t)]P[N{t)=1^^s,B{t)] ( 22 ) 

* j 

= + At)\s,N{t) = r]j,b{t)]P[N{t) = r]j\s,B{t)] (23) 

* j 

where we have used the Markov property P[iV(t + At) = rji^ b(t + At) |s, N(t) = r]j, B{t)] = P[A^(t + 
At) = rii,b(t + At)|s,iV(t) = rij,b{t)] to arrive at Equation (23). 


We now focus on the term P[iV(t + At) = rji, b{t + At)|s, iV(t) = rjj, b{t)] in Equation (23). This 
term is the state transition probability. Using the CTMP in Section we have 


where 


P[A(t + At) = rji, b{t + At)|s, N{t) = rjj^ b{t)] 

= bb{t+At),b{t) + lPl + bb{t+At),b{t)-lP2 + bb{t+At),b{t)P3 

Pi = K,Vj-iiR^Vj,R{M - b{t)) At 
P 2 = At 

P3 = rji^j 


(24) 


(25) 

(26) 
(27) 


where is the i?-th element of rjj, i.e. there are rjj^R signalling molecules in the receiver voxel, 
and 


Vi^j = 


where 


dij At if i ^ j 

1 — (Xrij^ji{M — b{t)) — jjib{t) — djj) At if i = j 


i^j 


(28) 


By substituting Equation (25) into Equation (23), we have 


where 


p[b{t + At)|s,i3(t)] — Sb{t+At),b{t) + lQl + <^6(t+Ai),b(i)-lQ2 + db{t+At) ,b{t)Q 3 


* 3 


(29) 


(30) 


( 31 ) 
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We will now determine Qi, Q 2 and Q 3 . 
For Qi, we have 


Qi = - b{t)) At P[A^(t) = rij\s,B{t)] 

i 3 

= A(M - h{t)) At YJYJ^v^,VJ-^RVj,RP[Nit) = r]j\s,B{t)] 

i 3 

= A(M - b{t)) At 2 = r]j\s, B{t)] 

3 s.t. 

= \{M-b{t)) At Y^r]j^RP[N{t) =r]j\s,B{t)] 

3 

= A(M — b{t)) At E[n/{(t)|s, B{t)\ 


(32) 

(33) 

(34) 


Note that in Equation (32), the sum is over all states 77 * with at least one signalling molecule in 


the receiver voxel, i.e. rjj^R ^ 1. Since the summand in Equation (32) is zero if r]j^R = 0, we get the 


same result if we are to sum over all possible states, that is why Equation (33) holds. 
For Q 2 , we have 


^ 3 


fib{t) AtY^YjK,rij+tR = Vj\-s,l3{t)] 

^ 3 

(35) 

fib{t) At2P[Ar(t) = 7]j\s,B{t)] 

3 

(36) 

fib{t) At 

(37) 


Note that Equation (36) follows from Equation (35) because for every r/j, there is a unique rji 
such that rji = r]j + Ijj holds. 
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For Q 3 , we have 


Q 3 = r]j\s,B{t)] 

i j 

=x;i;wj A*)p[A'(()=%k.B(()]+ 

i jj^i 

y~|(l — Xrij^ii{M — b{t)) At — iib{t) At — djj At)P[A^(t) = rij\s, B{t)] 

j 

= ^(1 - Xr]j^R{M - b{t)) At - iib{t) At)P[N{t) = r]j\s,B{t)] + 

j 

= Vj\s,l3{t)] -YjdjjP[N{t) = r]j\s,B{t)]) At 

i j^i j 

= (1 — A(M — 6(t))E[77,^(t)|s, B{t)\ At — ^ib{t) At) + 

(X;i;<i«p[A'(*) - %is.BW] -x;x;<i«p[iv(«) = At 

i j i^j 

'-V-' 

=0 

= (1 — A(M — 6(t))E[n/{(t)|s, B{t)] At — ^ib{t) At) (38) 


Having obtained Qi, Q 2 and Qs, we arrive at: 


P[ 6 (t + At)|s,H(t)] =5b{t+M),b{t)+iX{M -bit)) At E[nFt(t)\s,B{t)] + 

bb{t+/^t),b{t)-it‘d>{t) At + 

<5b(i+At),fe(t)(l - A(M - b{t))E[nR{t)\s, B{t)] At - ^b{t) At) 


Note that Equation (39) is the same as Equation (13) in the main text. 


(39) 


B Proof of Equation (14) 


From Equation (IT^, we have: 


dL,{t) ^ log(P[&(t + At)| 5 ,H(t)]) _ 

dt At^o At At^O 


log{P[b{t + At)\B{t)]) 
At 


(40) 


Note that the second term on the RHS is independent of transmission symbol s, we will focus on 
the hrst term. 


Note that P[b{t + At)\s,B{t)], which is given in Equation (39), is a sum three terms with 
multipliers 56(t+At),6(i)+i, db{t+At),b{t)-^ db{t+At),b{t)- Since these multipliers are mutually exclusive, 
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we have: 


log (P[&(t + At)\s, B{t)]) =Sb{t+At),b{t)+i log - Kt)) At E[n/j(t)|s, B{t)]) + 

Sb{t+At),b{t)-i log + 

Sb{t+At),b{t) log ((1 - A(M - 6(t))E[njj(t)|l, B{t)] At - ^b{t) At)) (41) 
^5b{t+At),b{t)+i log (E[n/j(t)|s, B{t)]) - 
<5b(t+At),6(t)A(M - 6 (t))E[nR(t)|s, B{t)] At+ 

Pit) (42) 


where we have used the approximation log(l + a At) ss a At in the last term of Equation (41) and 
have collected all terms that do not depend on s in Pit). 


By substituting Equation (42) into Equation (40), and taking limit At 0, we have 

dLsit) 


dt 


- iiSo ‘°e (EK(«)k,B(«)]) - 

- b(t)) (E[n„(i)|s,B(i)]]) + P) 

dUit) 


dt 


log(EK(t)|s,S(t)]) - A(M - bit)) (EK(t)|s,S(t)]) + Lit) 


(43) 

(44) 


where all terms that are independent of s have been collected in Lit). Note that Lit) contains some 
terms that diverges but this is not an issue because for demodulation it is their relative difference 
Lsi(f) — Lg^it) (for any two symbols si and S 2 ) that matters. 


Note also that we have used the following reasonings to arrive at Equation (44) from Equation 


(43): 


1. The term liiUAt^o is an impulse whenever a receptor changes from the unbound 

to the bound state. This is precisely 


2. The term 6 b{t+At),b{t) 1® o^ly when the number of bound receptors changes and the number 
of such changes is hnite. In other words, Sb{t+At),b{t) = ^ with probability one. This allows us 
to drop db(^t+At),b{t)- 


Finally, note that Equation (44) is the same as Equation (14). 


C An example end-to-end model 

For this example, we assume the transmission medium consists of 3 voxels as illustrated in Figure ??. 
The transmitter and receiver are assumed to be located in, respectively. Voxels 1 and 3. We assume 
reflecting boundary condition which means the signalling molecules cannot leave the medium. 
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This Appendix presents an example end-to-end model. For this end-to-end model, the trans¬ 
mitter is assumed to send Symbol 0 which means it uses the set of chemical reactions corresponding 
to this symbol. We therefore view a transmitter as a set of chemical reactions located within the 
transmitter voxel. It is still an open problem what chemical reactions are good for communication 
performance. The example being used here is not meant to promote the use of a particular set of 
chemical reactions but our purpose is to show how a set of chemical reactions can be modelled by 
a CTMP. 

For this example, we assume that the production of the signalling molecules S requires two 
intermediate chemical species F and G, which are produced by RNAi and RNA 2 . There are four 
reactions and they are assumed to take place within the transmitter voxel only. The four chemical 
reactions are: 


RNAi RNAi + F 

(45) 

RNA 2 RNA 2 + G 

(46) 

F S 

(47) 

S + G-^0 

(48) 


Reaction (45) says that the molecules of F are produced at a mean rate of ki. Similarly, according 


to Reaction (46), G is produced at a mean rate of ^ 2 - Reaction (47) says that F is converted to S at 


a mean rate equals to times the number of F molecules in the transmitter voxel. Reaction (48) 
says that S and G can react to produce a molecule 0 that we are not interested to keep track of in 


the mathematical model. If an S (or a G) molecule takes part in Reaction (48), we can consider 


this S (G) molecule has left the system permanently after the reaction. The rate of Reaction (48) 
is ^4 times the number of G molecules and the number of signalling molecules S in the transmitter 
voxel. 

We assume that the chemical species RNAi, RNA 2 , F and G are found in the transmitter voxel 
only, and they cannot leave the transmitter voxel. This means that we do not need to consider the 
diffusion of these chemical species. The only diffusible chemical species in the entire system is the 
signalling molecule S. We also assume that there is only one of each RNAi and RNA 2 and their 
counts remain constant. 

In order to dehne the state of the system, we make the following dehnitions: ni{t) is the number 
of signalling molecules in Voxel i at time t, npit) and ncit) are respectively the number of F and G 
molecules at time t, nA(t) is the cumulative number of molecules that have left the system at time 
t and b{t) is the number of complexes (or bound receptors) at time t. Since a receptor can either 
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be unbound or in a complex, the number of unbound receptors at time t is M — b(t)] therefore, 
the mathematical model only has to keep track of either the number of unbound receptors or the 
number of complexes, and we have chosen to keep track of the latter. The state of the system 
is completely specihed by these seven molecular counts: ni{t),n 2 (t),n 3 (t),nF(t),nG(t),nA{t) and 
b{t). All the molecular counts should be non-negative integers (i.e. belonging to the set Z^q) and 
a further restriction is that 0 < b(t) < M or we write b(t) e Z[o,m]- 
We dehne the vector N(t) as 


N{t) 


ni(t) n 2 (t) n^it) npit) ncit) UAit) 


T 


(49) 


where the superscript ^ is used to denote matrix transpose. 

Based on the dehnition of N{t), the state of the system is the tuple {N(t), b(t)) and a valid state 
must be an element of the set S = Z|g x Zp^M]- The state of the system changes when a reaction 
or diffusion event occurs. Our modelling assumptions mean that reactions can only take place in 


the transmitter or the receiver voxels. The reactions in the transmitter voxel are (45) —(48). The 
reactions taking place in the receiver voxel are 0 and ([^. The only diffusible chemical species in 
this system is the signalling molecule S. Within an inhnitesimal time At, at most one diffusion or 
reaction event can occur. Therefore, the dynamics of the system can be specihed by the transition 
probability from state {N(t),b{t)) to {N(t + At),b(t + At)). We will now specify these transition 
probabilities and we begin with the transmitter. 


Four possible reaction events (45)-(48) can take place in the transmitter voxel. An occurrence 


of Reaction (45) increases the number of F molecules in the transmitter voxel by 1 and this occurs 
at a mean rate of ki. By dehning Ij to be the standard basis vector with a ‘1’ at the Ath position. 


we can write the state transition probability due to Reaction (45) as: 


P[N{t + At) = N{t) + l4,b{t + At) = b{t)\N{t),b{t)] = ki At (50) 


Note that we have used I 4 because nF{t) is increased by 1 if Reaction (45) occurs and npit) is the 


fourth element of N{t) in the dehnition of N{t) in (49). The right-hand side (RHS) of Equation (50) 


is the transition probability that Reaction (45) occurs in {t,t + At), which is given by the reaction 
rate ki times At. 

We can write the transition probabilities due to Reactions (46) — ([4^ as: 


P[N{t + At) = N{t) + I 5 , b{t + At) = b{t)\N{t), b{t)] = fca At (51) 

P[A^(t -I- At) = N{t) — I 4 -I- li, b{t + At) = b{t)\N{t), b{t)] = ks npit) At (52) 

P[A^(f -I- At) = N{t) — I 5 — li -h 216, b{t + At) = b{t)\N{t), b{t)] = k^ n^t) ni{t)At (53) 


30 















The rationale behind Equation (51) is similar to that of (50). Equation (52) models Reaction (47). 


If Reaction (47) occurs, an F molecule is converted to an S molecule, so the number of F molecules 
(which is the fourth element of N{t)) is decreased by 1 and the number of signalling molecule 
in the transmitter voxel ni{t) (which is the hrst element of N(t)) is increased by 1; this change in 
the number of molecules as a result of Reaction ([4^ can be written as N{t + At) = N{t) — 1^ + li 


in (52). Equation (53) models Reaction (48). When Reaction (48) occurs, a G and an S molecule 


in the transmitter are consumed, hence —IL 5 — li in (53). We are not interested to keep track of 
the molecules as a result of this reaction, we consider these two molecules have left the system 
permanently and add ‘2’ to which is at the sixth position of N{t). The letter ‘A’ here comes 

from ’absorbing’ because once a molecule is added to nA(t), it will not leave. Note that the RHSs 


of (50)-(53) show the transition probabilities and they are of the form of the transition rate times 
At. 


The state of the system can also be changed by signalling molecules diffusing from one voxel 
to another. For this example, there are four possible diffusion events, which take place when a 
signalling molecule diffuses from a voxel to its neighbouring voxel. The four diffusion events are: 
from Voxel 1 to Voxel 2, from Voxel 2 to Voxel 1, from Voxel 2 to Voxel 3, and from Voxel 3 to 
Voxel 2. The transition probabilities of these four events are: 


P[V(t + At) = N(t) — li + IL 2 , b(t + At) = b(t)\N(t), b{t)] = d ni(t) At (54) 

P[V(t + At) = N{t) + li — I 2 , b{t + At) = b{t)\N{t), b{t)] = d n2{t) At (55) 

P[V(t + At) = N{t) — I2 + I3, b{t + At) = b{t)\N{t), b{t)] = d n2{t) At (56) 

P[V(t + At) = N{t) + I2 — I3, b{t + At) = b{t)\N{t), b{t)] = d n^{t) At (57) 


Equation (54) is the probability that a signalling molecules diffuses from Voxel 1 to Voxel 2. The 
occurrence of this event means the number of signalling molecules in Voxel 1 (= ni{t), which is the 
hrst element of N(t)) is decreased by 1 while the number of signalling molecules in Voxel 2 (= n 2 (t), 
which is the second element of N{t)) is increased by 1. The probability of this occurring is d At. 
The explanation for the other three transition probabilities are similar. 

The last category of state transitions occurs when a receptor is bound or unbound according to 
chemical reactions ([^ and (|^. The state transition probabilities are: 


P[N{t + At) = N{t) - I 3 , b{t + At) = b{t) + l|A^(t), b{t)] = A n^it) (M - b{t)) At (58) 
P[N{t + At) = N{t) + I 3 , b{t + At) = b{t) - l|A^(t), b{t)] = p b{t) At (59) 


Equation (58) is the transition probability for receptor binding or the formation of new complex. 


31 
















This event occurs when a signalling molecule in the receiver voxel reacts with a unbound receptor 
to form a complex. As a result of this reaction, the number of signalling molecules in the receiver 
voxel (which is Voxel 3 in this example) is decreased by 1 and the number of complexes b{t) is 
increased by 1. The rate of this event is proportional to the product of the number of signalling 


molecules in the receiver voxel n^it) and the number of unbound receptors {M — b(t)). Equation (59) 
is the transition probability for a receptor to unbind. The unbinding reaction causes the number of 
signalling molecules in the receiver voxel n^{t) to increase by 1 while the number of complexes b{t) 
to decrease by 1. The rate of this reaction is proportional to number of complexes b{t). 


Equations (50) to (59) give the transition probabilities of the possible events that can occur 
when the state of the system is {N{t),b(t)). It is possible that no transitions occurs in the time 
interval (t,t + At), the probability of this occurring is given by the complementary to that of an 


event occurring, that is, one minus the sum of the RHSs of Equations (50) to (59). This completes 
the model. 


D Difference between ^[nR{t)\s^ B{t)] and as{t) 


In Section]^ we propose to replace E[?7,jj(t)|s, i3(t)] by as{t) = E[?7,j^(t)|s]. The difference between 
these two quantities is hard to compute because ligand-receptor binding is a non-linear process. In 
this Appendix, we will provide some justihcation of this replacement by considering the hltering 
problem of LTI systems. Consider the following continuous-time LTI system in state space form 

dx{t) 


dt 


=Ax{t) + Bus{t) + w{t) 


y(t) =Cx{t) + Dus{t) -I- v{t) 


( 60 ) 

( 61 ) 


where (1) x{t), Ugit) and y{t) are, respectively, the state, input and output vectors; (2) w{t) is 
the state noise vector and v{t) is the measurement noise vector, where both vectors are zero-mean 
Gaussian white noise; and (3) (A, B, C, D) is a state space quadruple. We assume the dimensions of 
the vectors and matrices are compatible. The subscript s in Us{t) plays the same role as the Symbol 
s in the main text. We can use different s to choose different input Us{t). We can determine the 


mean state vector E[x(t)] by taking expectation on both sides of Equation (60) to obtain: 

(iE[a;(t)] 


dt 


=AE[a:(t)] -I- Bus{t) 


(62) 


Note that E[x(t)] plays the same role as as(t) = E[ni?(t)|s]. 
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The filtering problem for a LTI system is to estimate the state vector x{t) from the continuous 
history of the output y(t). A method to realise filtering is to use an observer [6T] : 

dx(t) 


dt 


=Ax(t) + Bu,{t) + K(y(t) - p(t)) 


y(t) ^Cx(t) + Du^(t) 


(63) 

(64) 


where K is the observer gain matrix. The vector x(t) is the estimated state vector from the 
past history of the output y(t). The expectation of x(t), i.e. E[a;(f)], plays the same role as 
E[nR(t)js,jl3(t)]. 


We are interested to study the difference e{t) = E[a;(f)] — E[x(f)]. By using Equations (60), 
ij), ( [6^ and (64), it can be shown that e{t) obeys the differential equation 

de{t) 


dt 


= {A + KC)e{t) 


(65) 


This means that if the matrix A + KC is asymptotically stable, then the difference e(t) tends to 
zero for sufficiently large t. This result suggests that the difference between E[nij(t)|s, i3(t)] and 
as{t) can be small if the filtering error is stable. 
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Figure 1: An overview of the system considered in this paper. 



Figure 2: A model of molecular communication network. The volume is divided into voxels. The 
indices of the voxels are given in the top right hand corner. Unhlled circles are signalling molecules. 
Filled circles are receptors. The dashed lines show the boundary of the medium. 
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Figure 3: This figure explains the meaning of the function U{t), which is the total number of 
unbound-to-bound transitions for all receptors. 



Figure 4: The demodulator structure. 
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Figure 5: The figure compares (7s(t) = E[ni^(t)|s] and E[nij(t)|s, i3(t)] for one realisation of B{t). 
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Figure 6: The SERs of the optimal and sub-optimal filters for different decoding time t. 
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Symbol 0 is sent 



Figure 7: The output of the modulators Z^it) (thin line) and Zi{t) (thick line) for Symbol 0. 


Symbol 1 is sent 



Figure 8: The output of the modulators Zo(t) (thin line) and Zi(t) (thick line) for Symbol 1. 
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Symbol 0 is sent 



Figure 9: The mean output of the modulators ^o(^) (thin line) and Zi{t) (thick line) for Symbol 0. 



Figure 10: The SER for Symbols 0 and 1. 
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Figure 11: The SER for Symbols 0 and 1 for varying number of receptors. 



Figure 12: The SER for Symbols 0, 1 and 2 for varying number of receptors. 
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Symbol 0 is sent 



Figure 13: The output of the modulators ^o(^) (thin line) and Zi{t) (thick line) for Symbol 0. The 
mean number of signalling molecules at the receiver voxel for both symbols is similar. 



Figure 14: The average SER versus number of receptors. The ISI case. 
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